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I. INTRODUCTION 

The enormous advances in parallel computing made 
during the past few years, together with theoretical ad- 
vances in the formulation of lattice gauge theories with 
fermions, have allowed lattice theorists to abandon the 
quenched approximation that dominated lattice QCD 
simulations for such a long time in favour of simulations 
using dynamical light quarks. This important step has 
allowed a significant reduction in systematic errors by re- 
moving the large and uncontrolled errors inherent in the 
quenched approximation. 

The Fermilab Lattice, MILC and HPQCD collabora- 
tions have an ambitious program which to date has made 
several high-precision predictions from unquenched lat- 
tice QCD simulations [l], 0], including accurate deter- 
minations of the strong coupling constant Ug @, the 
light and strange quark masses [j] , and the leptonic and 
semileptonic decays of the D meson [Ej. To do this, we 
rely on the Symanzik-improved staggered-quark formal- 
ism [6,] , specifically the use of the asqtad [3] action. While 
this approach requires the use of the fourth root of the 
staggered quark action determinant, all of the available 
evidence to date is consistent with the conclusion that 
the resulting theory is in the same universality class as 
continuum QCD, as long as the chiral limit is taken after 
the continuum limit Q. 

Recent studies of the heavy-quark potential in full 
QCD [9| have shown an apparent increase in scaling vi- 
olations compared to the quenched approximation, con- 
trary to expectations. A possible reason for this would be 
that these scaling violations arise from the mismatch be- 
tween the inclusion of sea quark effects in the simulation 
and the omission of sea quark effects in the improvement 
coefficients in the action, which would appear to spoil 
the 0{a^) improvement at the level of 0{asN fo?). A 
systematic study of 0{asa?') effects is generally beyond 
the scope of the current perturbative improvement pro- 
gramme. Nevertheless, it is important to bring up-to- 
date the calculation by Liischer and Weisz jl^l and by 



Snippe [Tlj of the radiative correction to the 0{a?) tree- 
level Symanzik-improved gluon action to include the ef- 
fects of dynamical quarks. This is important also because 
the Liischer-Weisz improvement is currently included in 
many unquenched simulations 0] • Since the lattice spac- 
ing scale is set by measurement of the heavy-quark po- 
tential, there will be an induced 0{asN fo?') artifact by 
omitting the corrections due to unquenching. While such 
errors are generally smaller than other systematic errors 
in current state-of-the art studies, it is simple to remove 
them, using the result of the perturbative matching cal- 
culations done here, and this may prove advantageous in 
careful studies of different scale setting procedures. 

In this paper, we present the determination of the 
lowest-order perturbative contributions from quark loops 
to the Symanzik improvement coefRcients of the Liischer- 
Weisz glue action. Including these contributions in fu- 
ture simulations, as well as accounting for their influ- 
ence in the analysis of existing results, should help to 
eradicate the last remaining vestiges of the quenched ap- 
proximation and any associated systematic errors from 
unquenched lattice results. Some of this work has been 
reported in preliminary form in [l^ . 



II. CONCEPTS AND METHODS 

First, let us briefly explain the ingredients of our cal- 
culation. 



A. On-shell improvement 

The original Symanzik improvement programme fl^, 
[T^ aims to remove the discretisation artifacts from the 
correlation functions of the lattice theory. For gauge the- 
ories, this has proven difficult to implement, since the 
correlation functions themselves are not gauge invariant. 
A way out of this difficulty is offered by the method of 
on-shell improvement introduced by Liischer and Weisz 
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[isl . [l^ which aims to improve only gauge invariant spec- 
tral quantities. 

The Liischer-Weisz action is given by [H, [13] 



(1 - P^,) + 2ci V (1 - R,,) (1) 



where P, i? and T are the plaquette, rectangle and 
"twisted" parallelogram loops, respectively. 

The requirement of obtaining the Yang-Mills action in 
the continuum limit imposes the constraint 



Co + 8ci + 8c2 = 1 



(2) 



which can be used to determine cq in terms of the other 
two coefficients. This leaves us with ci and C2 as un- 
known coefficients which need to be determined in order 
to eliminate the C(a^) lattice artifacts. 

If we have two independent quantities Qi and Q2 which 
can be expanded in powers of (/ia), where fi is some en- 
ergy scale, as 

= + + O (i^iay) (3) 

and which receive corrections 

Ai„,pQ, = d,^cj{fiaf + O {{^laf) (4) 

from the improvement operators, then the 0{a^) match- 
ing condition reads 



■'W^ 



(5) 



Since this equation is linear, we can decompose the Wi 
into a gluonic and a fermionic part as Wi = it;f'"° + 
Nfwf^^^''^ and obtain the same decomposition for the c^; 
thus, especi ally we do not need to repeat the quenched 
calculation [lH in order to obtain the 0{Nf) con- 
tributions (however, doing so provides a useful check on 
the correctness of our methods, which we have performed 
successfully). At higher orders in perturbation theory, 
the dij and Wi will become functions of the Ci in lower 
orders. 

At the tree-level, the fermions contribute nothing to 
gluonic observables, and hence the tree-level coefficients 
remain unchanged compared to the quenched case [Toj : 



1 

l"2 



C2 = 0. 



B. Lattice perturbation theory 



(6) 



Lattice field theory is usually employed as a non- 
perturbative regularisation; for the calculations we need 



to perform, however, we need a perturbative expansion 
of Lattice QCD. 

In lattice perturbation theory, the link variables are 
expressed in terms of the gauge field as 



U^{x) ^ exp {gaA^ (x + Ifi)) 



(7) 



which, when expanded in powers of g, leads to a per- 
turbative expansion of the lattice action, from which the 
perturbative vertex functions can be derived. 

The gauge field is Lie algebra- valued, and can be 
decomposed as 



(8) 



with the being anti-hermitian generators of SU(Af), 
where = 3 in the case of QCD. 

As in any perturbative formulation of a gauge theory, 
gauge fixing and ghost terms appear in the Fadeev-Popov 
Lagrangian; here we will not have to concern ourselves 
with these, since for the purpose of determining the un- 
quenching effects at one loop we only need to consider 
quark loops. An additional term, which we also do not 
need to consider here, arises from the Haar measure on 
the gauge group. 

The loop integrals of continuum perturbation theory 
are replaced by finite sums over the points of the recip- 
rocal lattice in lattice perturbation theory, or integrals 
over the Brillouin zone where the lattice has infinite spa- 
tial extent. 

To handle the complicated form of the vertices and 
propagators in lattice perturbation theory we employ a 
number of automation methods [3, [l^, llOj HH 113 ^^^^ 
are based on the seminal work of Liischer and Weisz [l^ . 
Three independent implementations by different authors 
have been used in this work to ensure against program- 
ming errors. 



C. Twisted boundary conditions 

We work on a four-dimcnsional Euclidean lattice of 
length La in the x and y directions and lengths LzU, Lta 
in the z and t directions, respectively, where a is the 
lattice spacing and L,Lz,Lt are even integers. In the 
following, we will employ twisted boundary conditions 
(23j in much the same way as in [l^. El . The twisted 
boundary conditions we use for gluons and quarks are 
applied to the (x, y) directions and are given by {v = x, y) 



Ut,{x + Lv) = n^Uf,{x)Q^- 



(9) 
(10) 



where the quark field ^scix) becomes a matrix in smell- 
colour space [13] by the introduction of a "smell" group 
SU(A^s) with Ns = N in addition to the colour group 
SU(A'^). We apply periodic boundary conditions in the 
(z, t) directions. 
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These boundary conditions lead to a change in the 
Fourier expansion of the fields which now reads 



^ p 



(12) 



In the twisted (x, y) directions the sums are over 



NL NL 



^ = {x,v) , (13) 



where m = and in the untwisted (z, t) directions the 



sums are over 
27r 



Pu 



2 



{z,t). (14) 



The modes with (rix — ny = mod N) are omitted 
from the sum in the case of the gluons since the gauge 
group is non-abelian, and this is signified by the prime on 
the summation symbol in Eqn. (Ilip . In particular, this 
removes the zero mode from the gluon spectrum and so 
the mass-scale m defined above acts as a gauge-invariant 
infra-red regulator. The matrices Tp are given by (up to 
an arbitrary phase, which may be chosen for convenience) 



(15) 



The momentum sums for quark loops need to be divided 
by N to remove the redundant smell factor. 

The twisted theory can be viewed as a two-dimensional 
field theory in the infinite {z,t) plane with the modes in 
the twisted directions being considered in the spirit of 
Kaluza-Klein modes. Denoting n — {nx,ny), the stable 
particles in the {z,t) continuum limit of this effective 
theory are called the A mesons (n ~ (1, 0) or n = (0, 1)) 
with mass m and the B mesons (n = (1, 1)) with mass 
V2m [ll|. 



D. Small- mass expansion 

To extract the 0{a'^) lattice artifacts, we first expand 
some observable quantity Q in powers of ma at fixed mqa: 

Q(ma,mqa) — a'^^niqa) + a'^\mqa){ma)'^ + 

O {{jnaY, [maY log(ma)) (16) 

where the coefficients in the expansion are all functions 
of mqa. There is no term at O [{ma)^ log(raa)) since 
the gluon action is improved at tree-level to O(a^) [ll[. 
Since we wish to extrapolate to the chiral limit it might 
be thought that we can set m^a = straight away to 
achieve this end. However, the correct chiral limit is 
mqa 0, ma 0, mq/m > C, where m = as be- 
fore and C is a constant determined by the requirement 



that the appropriate Wick rotation can be performed in 
order to evaluate the Feynman integrals. If the inequal- 
ity is violated this results in a pinch singularity. It is 
physically sensible that the correct limit is L — > oo be- 
fore mq — > since this divorces the two infra-red scales 
and avoids complication. This does, however, require us 
to consider the double expansion in m^a, ma and carry 
out the extrapolation to mqa = for the coefficients in 
Eqn. p6)) . We return to this issue in the next section 
when we discuss choice of integration contours. 

To extrapolate to the chiral limit, m^a — > 0, we will 
fit the coefficients in the expansion for Q in ma to their 
most general expansion in mqa for small m^a. 

For a^\mqa) we have 



AQ) 



(mqa) = b'^^fj log(m,a) + a'^^^ 



(Q) 



(17) 



Since we expect a well-defined continuum limit, 
a^^\mqa) cannot contain any negative powers of m^a 
but, depending on the quantity Q, it may contain loga- 
rithms; y^Q^Q is the continuum anomalous dimension as- 
sociated with Q which is determined by a continuum cal- 
culation. 

There can be no terms in [mqa)"^"" ,n > since these 
are obviously non-zero in the limit ma 0, and there is 
no counterterm in the gluon action that can compensate 
for a scaling violation of this kind. 

AQ) 



For 02 (wqa) we find 



AQ) 



AQ) 



(mqa) 



{mqaf 



AQ) 



(18) 



+ (02^2 + 4? log("iga)) {mqaf + O ((m,a)^ 



After multiplication by [ma)"^ the (mqa) ^ contribution 

gives rise to a continuum contribution to Q, and 02*^2 is 
calculable in continuum perturbation theory. There can 
be no term in (ruqa)^^ \og{mqa) since this would be a 
volume-dependent further contribution to the anomalous 
dimension of Q, and there can be no term in \og(mqa) 
since the action is tree-level O(a^) improved [13]. 

Depending on the choice of observable Q there may be 
additional constraints on the coefficients which appear in 
the expansions. We discuss these in the next section in 
the context of the particular observables with which we 
concern ourselves. 



III. CALCULATIONS AND RESULTS 

In this section we lay out the calculation of the un- 
quenching effects to order 0{asNfa^). The numbers and 
quantities given in the following are per quark flavour, 
and hence need to be multiplied by Nf throughout. 
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Lorentz invariance will be restored, although for L finite 
this will not be the case. Gauge invariance and the Ward 
Identity then ensure that, in this limit, Lorentz invari- 
ance implies that the gluon self-energy function 7rp^(k) 
satisfies 



(22) 



FIG. 1: The fermionic one- loop diagrams contributing to the 
A meson mass renormalisation as well as to the wavefunction 
renormalisation for A and B mesons 



A. The A meson mass 

The simplest spectral quantity that can be chosen 
within the framework of the twisted boundary conditions 
outlined above is the (renormalised) mass of the A me- 
son. In agreement with Eqn. (109) of [llj the one-loop 
correction the the A meson mass (for A mesons with pos- 
itive spin) is given by 



,(1) 



-^o(k) 



^«(fc) 



2m 



(0) 



(19) 



k={im^2 fi.mfi) 



where .^o(k) = 1 + ((ma)"*) is the residue of the pole 
of the tree-level gluon propagator at spatial momentum 
and TO^' is defined so that the momentum k is on- 



shell. We consider the dimensionless quantity rn'^' /m. 
The fermionic diagrams that contribute to this quantity 
are shown in figure [TJ 

The anomalous dimension of uia is zero and so using 
Eqn. p?]) we have 



,(m^,l) _ 







(20) 



Using physical arguments we can determine the be- 
haviour of other coefficients. From continuum calcula- 
tions we find 



a 



("IA,1) 

2,-2 







(21) 



This result follows from the fact that the fermion con- 
tribution at one loop order to niA is IR finite since the 
fermion has a non-zero mass and is 4D Lorentz invariant. 
Thus, a'^_^2 can be constructed only from 4D Lorentz in- 
variants of which we have only ea ■ k^ and k^, where ea 
is the j4-meson polarization vector. However, gauge in- 
variance implies • k^ = and the on-shell condition 
gives k^ = and so, there being no non-zero Lorentz 
invariant, we deduce the result. 

A much less obvious deduction is that a^^'^\mqa) — 
0, which together with eqns. (|17I20[) implies that 

flg™'*'^'' = 0. A necessary ingredient to derive this result 
is the fact that the one-loop fermion contribution to tua 
is IR finite in the limit m — > (L — > oo) since the fermion 
mass, m„, is non-zero. We thus expect that in this limit 



From [id . lll| the one- loop contribution to iriA is propor- 
tional to e'^e''7r^^(k). In the limit L ^ oo we are able to 
use Eqn. and we find that this contribution is zero 
by gauge invariance and the on-shell condition fc^ = 0. In 
contrast, the contribution to ag™'^'* from internal gluon 
loops is not zero by this argument; it is indeed calculated 
in \Td, Tl'l . The reason is that the one- loop gluon contri- 
bution is not IR finite in the to ^ 0, L — > oo limit since 
the IR-regulating mass is m; the internal gluon "feels" 
the finite boundary of the lattice in the x, y direction no 
matter how large L is. Consequently, we cannot expect 
to use restoration of Lorentz invariance to limit the form 
that the purely gluonic 7r^i/(k) takes, and so no deduction 
concerning this coefficient can be made. 

An alternative explanation for why Oq'"'*'^'' — also re- 
lies on the restoration of Lorentz invariance in the m — s- 
limit. In this limit, the action is isotropic with metric ten- 
sor g^i, = diag(l, 1, 1, 1). However, the twisted boundary 
conditions break Lorentz invariance and single out the 
twisted X, y directions, and so we must expect that ra- 
diative corrections will renormalise g^^i, in a way that can 
break Lorentz symmetry. The mass shell condition for 
the A-mesou is then 



g'l^^k^'k" = 0, fc^ = {ipo, 0, to, fcg) 



(23) 



where g^ is the renormalised metric tensor. This is rein- 
terpreted as a renormalization of the j4-meson mass to 
with 



511 



(24) 



This can also be interpreted as an anisotropy renormal- 
ization. Since the one-loop fermion contribution is IR fi- 
nite and Lorentz symmetry is restored in the limit to — > 
we then have that gii{m = 0) = gu = 1 and wla is 
not renormalised. This is not the case for the one-loop 
gluon contribution, which is not IR finite, and so the as- 
sumption that Lorentz symmetry is restored as to —> is 
incorrect. 

For the kinematics used here this means that in the 
limit TO ^ 0, L oo then tth vanishes and hence from 
Eqn. (|19p so does rrS^ Im. This expectation is accurately 
verified by our calculation: the extrapolation of to^-'/to 
to TO = indeed gives zero (cf. figure [31). 

In the chiral limit ruq — > 0, the term Wi that appears 

on the right-hand side of Eqn. ([5]) is aj*^"* , and it is this 
limit and this coefficient that we will concern ourselves 
with hereafter. 



5 



The O {ctsirna)^^ contribution from improvement of 
the action is given by [ll| 



,(1) 



-(4'^ - 4'^)(ma)2 ^ Q [{maf) (25) 



leading to the improvement condition 



(1) _ (1) _ 
Cl C2 — '^2.0 



B. The three-point coupHng 



(26) 




An effective couphng constant A for an AAB meson 
vertex is defined as 



A - .goV^(k)Z(p)Z(q)e,ri'2^^-(fc,p,g) (27) 

where we have factored out a twist factor of 
Tr([rfc, Fpjrg) from both sides, and the momenta and 
polarisations of the incoming particles are 

fc = (i£;(k),k); k=(0,m,«r) 

p = (-i£:(p),p); p=(m,0,zr) 

g=(0, q); q = (— m, — m, — 2ir) 
e=(0, 1,-1,0) (28) 

Here r > is defined such that -E(q) = 0. This cou- 
pling is a spectral quantity since it can be related to the 
scattering amplitude of A mesons We expand Eqn. 
([?f)) perturbatively to one-loop order and find in agree- 
ment with Eqn. (137) of fuj 



aw 

m 



1 2- 

24"^ , 



r(i) 



/co=i-B(k 



(29) 



The fermionic diagrams contributing to the irreducible 
three-point function V^^^ are shown in figure [21 Using 
Eqn. ()17|) and the known anomalous dimension of the 
coupling constant we have 



"0,0 



f „2 



37r2 



(30) 



Unlike the argument for Cj"!'^'^'' 
calculation gives 



above, a continuum 



1207r2 



(31) 



In this case there are non-zero Lorentz invariants for the 
three-point function such as -k^ etc. and so we expect 
this coefficient to be non-zero. 



The improvement contribution to A is 11 1 
Ai 



m 



4(9ci'^ - 74'^)(mo)2 + O ((ma)^) (32) 



FIG. 2: The fermionic one- loop diagrams contributing to the 
three-point function. 



leading to the improvement condition 

4(9c«-7c«) = -4^oi) (33) 



Tests of our calculation are that the fit for a, 



0.0 



must 



give the correct anomalous dimension stated in Eqn. 
pop , and that our fits reproduce the continuum result 

a2^^2 — — g^/1207r. Both are accurately verified (cf. fig- 
ures m and [5]) . 



C. Choice of integration contours 

The external lines of the diagrams are on their respec- 
tive mass shells but with complex three-momentum k in 
which the third component, fca, has been continued to an 
imaginary value parametrised by the variable r as shown 
in Eqn. ([28]); in the Euclidean formulation fep is also 
imaginary. In evaluating the loop integrals that are not 
pure tadpoles, care must be taken to ensure that the am- 
plitudes calculated are the correct analytic continuations 
in r from the Minkowski space on-shell amplitudes de- 
fined with real three-momenta to the ones in Eqn. ([^5)) . 

The situation is complicated by the presence of two 
mass scales m, m^. The integrals are evaluated after per- 
forming a Wick rotation in ko, taking care to avoid con- 
tour crossing of any poles that move as r is continued 
from r = to r = m/V2. We find that mq/m must 
be chosen larger than a minimum value, dependent on 
the graph being considered, to avoid any contour being 
pinched. The outcome is that after the Wick rotation in 
/co, the (Euclidean) integration contour for either or, 
in one case, k^ must be shifted by an imaginary constant. 

,(1) 



We find it sufficient that for the calculation of 



and Zyi(fc) we impose > m/2 and for F*-^) and Zsiq) 
that niq > raj\f2. 



6 



7W 



0.15 0.0036752(7) 

0.2 0.003701(1) 

0.3 0.003730711(1) 

0.4 0.00372996(4) 

0.5 0.003696507(2) 

0.6 0.0036328671(4) 

0.7 0.0035435429(3) 

0.8 0.0034337690(4) 

0.9 0.0033087971(2) 

1.0 0.0031734700(3) 

1.2 0.0028882(3) 



0.07456(1) 

0.0648161(5) 

0.051090(2) 

0.0414332(1) 

0.03408572(2) 

0.028272563(9) 

0.023575013(4) 

0.019733513(1) 



-0.178(6) 

-0.1617(4) 

-0.1498(9) 

-0.14498(6) 

-0.140933(7) 

-0.136776(2) 

-0.132222(1) 

-0.12722(3) 



0.009976(2) -0.1044(1) 



TABLE I: The coefficients from the fits of mj^' and A*^' /m 
against ma. 



Triple -Gluon: g ^( arn^) 




D. Extracting the coefficients 

To extract the improvement coefficients from our dia- 
grammatic calculations, we compute the diagrams for a 
number of different values of both L and with Nf = 1, 
N = .i. At each value of m^, we then perform a fit in ma 
of the form given in Eqn. (|16p to extract the coefficients 
{iriqa), n — 0,2. The results of these fits are given 
in Table m 

To facilitate our fits, we make use of the prior phys- 



we have 



ical knowledge we have: In the case of 
_ Q ]-,g(2ayse of gauge invariance. 
Performing a fit of the form ((T7]) and ((TS]) , respectively, 
on these coefiicients, we get the required coefficients of 
the Oia^) lattice artifacts in the chiral limit to be 



a 



(mA,l) 



2,0 



0.00361(1) 
-0.140(1) 



(34) 
(35) 



FIG. 4: A plot of flg^'^' against niqa which shows the agree- 
ment between the numerical lattice results and the known 
anomalous dimension. 



Here, again, we have facilitated our fits by making 
use of our prior knowledge: For a^2^_^2^'^ vanishes, 

and for \^^\ we have two known continuum contribu- 



tions: 6, 



'0,0 



-I/Btt^ is the one- loop coefficient of the 



/3-function and a2^^l — —XjYlOit^ is the continuum co- 
efficient of the infrared divergence m? /m?^. 

Solving eqns. (|^ and ([55)1 for c|^\ our results can be 
summarised as 



-1 



-0.025218(4) -f 0.00486(13)iV/ (36) 
-0.004418(4) + 0.00126(13)iV/ (37) 



These coefficients are to be identified with the Wi of Eqn. 



where the quenched {Nf — 0) results are taken from [Tlj. 
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FIG. 3: A plot of the fermionic contributions to the one-loop 
A meson self-energy rrS^ jm against (maf' . The vanishing of 
rrij^ /m in the infinite- volume limit can be seen clearly. 



FIG. 5: A plot of 0^2^'^^ against m^a with the analytical con- 
tinuum result for the infrared divergence shown for compari- 
son. 
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IV. CONCLUSIONS 

Repeating the analysis of JJ] and using their notation 
we exp ress the radiatively corrected action of Eqn. ^ 
as 



S[U] = ^ft5,[C/] 



(38) 



Then 



f3i 



12n 



i=0 



,(1) 



- 487rci 



(1) 



2u, 



(1) 



UttPo (1) 



5w 



(39) 



The quenched radiative contributions have been analyzed 
in [l7[ and so we may write 



/3i = - 



f3o 



1 + 0.4805a^ - ( ^4^1 + 487rc^^| | a. 



^0 f 0.033a, - i|^4>. 



(40) 



where now all the one-loop coefficients cf'j contain only 
quark loop contributions. 

Plugging in the numbers obtained in this work we find 



/3i = 



P2 



/3o 



20u2 



[f + 0.4805as - 0.3637(f 4)7V/as] , 



(0.033a, - 0.009(l)iV/a,) 



(41) 



With Nf — 3 the shift from the quenched values is sur- 
prisingly large, and may have a significant impact; espe- 
cially, it may explain the increased scaling violations seen 
in some unquenched simulations. 
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